MuskingumCunge Subroutine

SUBROUTINE MuskingumCunge &
!
( dt, dx, Qin, Pin, Pout, Qlat, Plat, B0, alpha, slope, n, Qout, depth, width, k, x)

!arguments with intent in:
INTEGER (KIND = short), INTENT (IN) :: dt !!time increment (s)
REAL (KIND = float), INTENT (IN) :: dx !!reach length (m)
REAL (KIND = float), INTENT (IN) :: Qin !!input discharge at current time (m3/s)
REAL (KIND = float), INTENT (IN) :: Pin, Pout !!input and output at previous time (m3/s)
REAL (KIND = float), INTENT (IN) :: Qlat, Plat !!lateral flow at current and previous time (m3/s)
REAL (KIND = float), INTENT (IN) :: B0 !!bottom width, = 0 for  triangular  section (m)
REAL (KIND = float), INTENT (IN) :: alpha !!angle formed by dykes over a horizontal plane (deg)
REAL (KIND = float), INTENT (IN) :: slope !!bed slope (m/m)
REAL (KIND = float), INTENT (IN) :: n !!manning roughness coefficient (s m^(-1/3))
REAL (KIND = float), OPTIONAL, INTENT(IN) :: k !!flood wave travel time used to compute C1, C2, and C3 (s)
REAL (KIND = float), OPTIONAL, INTENT(IN) :: x !!used to compute C1, C2, and C3 (-) 


!arguments with intent out:
REAL (KIND = float), INTENT (OUT) :: Qout !!output discharge at current time
REAL (KIND = float), INTENT (OUT) :: depth !average water depth 
REAL (KIND = float), INTENT (OUT) :: width !wetted width

!local declarations

REAL (KIND = float) :: xx, kk !!used for computing C1, C2, and C3
REAL (KIND = float) :: cel !!kinematic celerity (m/s)
REAL (KIND = float) :: normal_depth !!normal flow depth (m)
REAL (KIND = float) :: Qref !!reference discharge at time t+deltat
REAL (KIND = float) :: B ! topwidth
REAL (KIND = float) :: C1, C2, C3, C4
REAL (KIND = float) :: storage !storage at time t + deltat
REAL (KIND = float) :: aveWetArea !average wetted area in the reach

!-------------------------------end of declaration-----------------------------
IF (PRESENT (x) .AND. PRESENT (k) ) THEN
	  C1 = (dt-2.*k*x) / (2.*k*(1.-x) + dt)
	  C2 = (dt+2.*k*x) / (2.*k*(1.-x) + dt)
	  C3 = (2.*k*(1.-x) - dt) / (2.*k*(1.-x) + dt)
      C4 = (2.*dt) / (2.*k*(1.-x) + dt)
      Qout = C1 * Qin + C2 * Pin + C3 * Pout + C4 * 0.5 * (Qlat + Plat)
      
      
      
     !Compute the storage
      storage = k * x * Qin + k * (1. - x) * Qout
 
    IF (storage > 0.) THEN
        !Compute the average wetted area
        aveWetArea = storage / dx
 
        !Compute average water depth
        depth = DepthFromArea (aveWetArea, B0, alpha)
        
        !Compute width
        width = TopWidth ( depth, B0, alpha )

    ELSE
         depth = 0.
         width = 0.
    END IF
    RETURN !no need  to go on
END IF


!compute reference discharge
Qref = (Qin + Pin + Pout) / 3.

IF (Qref == 0. .AND. Qlat > 0.) THEN
  Qref = Qlat
END IF

IF (Qref > 0.) THEN
  
    !compute normal flow depth
    normal_depth = NormalDepth ( Qref, B0, alpha, slope, n )

    !compute topwidth
    B = TopWidth ( normal_depth, B0, alpha )

    !compute celerity
    cel = Celerity ( normal_depth, B0, alpha, slope, n )

    !compute k and x
    kk = dx / cel
    xx = 0.5 * ( 1. - Qref / (B * cel * slope * dx) )
    
    !compute weigths
	  C1 = (dt-2.*kk*xx) / (2.*kk*(1.-xx)+dt)
	  C2 = (dt+2.*kk*xx) / (2.*kk*(1.-xx)+dt)
	  C3 = (2.*kk*(1.-xx)-dt) / (2.*kk*(1.-xx)+dt)
      C4 =  2.*dt / (2.*kk*(1.-xx)+dt)
    
    !compute Qout
    Qout = C1 * Qin + C2 * Pin + C3 * Pout + C4 * (Qlat+Plat) / 2.
    !Qout = C1 * Qin + C2 * Pin + C3 * Pout + C4 * Qlat
    
    !last check
    IF (Qout < 0.) then
          Qout = 0.
    END IF
    
    !compute the storage
    storage = kk *xx * Qin + kk * (1. - xx) * Qout
    
    !compute the stage
    IF (storage > 0.) THEN
        !Compute the average wetted area
        aveWetArea = storage / dx
 
        !Compute average water depth
        depth = DepthFromArea (aveWetArea, B0, alpha)
        
        !Compute width
        width = TopWidth ( depth, B0, alpha )
        
    ELSE
        depth = 0.
        width = 0.
    END IF


ELSE
  Qout = 0.
  depth = 0.
END IF      

RETURN
 
END SUBROUTINE MuskingumCunge